% This practical we will work with a single subject's data from a blocked
% finger tapping experiment (CTF data courtesy of Matt Brookes), and
% perform a time-frequency analysis in the beta band using a time-wise GLM
% in source space.   
 
% Data can be downloaded from: 
% www.fmrib.ox.ac.uk/~woolrich/ft_subject1_data.tar.gz

%%%%%%%%%%%%%%%%%%
%% SETUP THE MATLAB PATHS
% make sure that fieldtrip and spm are not in your matlab path

global OSLDIR;
    
% set this to where you have downloaded OSL and the practical data:
practical_dir='/home/mwoolrich/Desktop'; 
osldir=[practical_dir '/osl1.3.1'];    

%practical_dir='/Users/woolrich';
%osldir=[practical_dir '/homedir/matlab/osl1.3.1'];    

addpath(osldir);
osl_startup(osldir);

%%%%%%%%%%%%%%%%%%
%% INITIALISE GLOBAL SETTINGS FOR THIS ANALYSIS

% directory where the data is:
workingdir=[practical_dir '/ctf_fingertap_subject1_data']; % directory where the data is
%workingdir=[practical_dir '/homedir/matlab/osl_testdata_dir/ctf_fingertap_subject1_data']; % directory where the data is

cmd = ['mkdir ' workingdir]; if ~exist(workingdir, 'dir'), unix(cmd); end % make dir to put the results in

cd(workingdir);

% Set up the list of subjects and their structural scans for the analysis 
% Currently there is only 1 subject.
clear spm_files;
spm_files={[workingdir '/dsubject1.mat']};
         
structural_files = {[workingdir '/subject1_struct.nii']};      

cleanup_files=0; % flag to indicate that you want to clean up files that are no longer needed

%%%%%%%%%%%%%%%%%%%%
%% load the file in and take a look at the SPM object.

subnum = 1;             
D = spm_eeg_load(spm_files{subnum});

% look at the SPM object. Note that it is continuous data, with 232000 time
% points at 250Hz. We will epoch the data later.
D

%%%%%%%%%%%%%%%%%%%
%% check data in oslview

oslview(D);

%%%%%%%%%%%%%%%%%%%
%% setup osl

oat=[];
oat.source_recon.D_continuous=spm_files;
oat.source_recon.conditions={'Undefined'};
oat.source_recon.freq_range=[4 30]; % frequency range in Hz
oat.source_recon.time_range=[300,32*30];
%oat.source_recon.time_range=[300,12*30];
oat.source_recon.method='beamform';
oat.source_recon.gridstep=10; % in mm, using a lower resolution here than you would normally, for computational speed
oat.source_recon.mri=structural_files;
do_hmm=0;
if(do_hmm)
    oat.source_recon.hmm_num_states=10;
end;

oat = osl_check_oat(oat);

%%%%%%%%%%%%%%%%%%%%%%%
%% run beamformer

oat.to_do=[1 0 0 0];
oat.source_recon.pca_dim=150;
oat.source_recon.dirname=[workingdir '/subj1_results_beta_hmm' num2str(do_hmm) '_' oat.source_recon.method];

oat = osl_run_oat(oat);

%%%%%%%%%%%%%%%
%% Establish regressor for continuous time GLM. 
% This should be setup to correspond to the time window
% over which the source recon is carried out.

D=spm_eeg_load(spm_files{1});

% This should be setup to correspond to the same time window
% as the full time window for D

x=zeros(length(D.time),5);

block_length=30; %s
block_order=[5 5 5 5 5 5 5 5 5 5 4 3 2 1 2 3 1 4 3 4 1 3 2 1 4 4 2 1 3 3 4 1 4 3 1 2 1 2 3 4 3 4 1 2 3 4 1 2];

% [Left, Right, Rest, Both, Rest_at_start]
% [  1     2      4     8     16 ]
% figure;plot(D.time,squeeze(D(1,:,:)))
% emacs ~/vols_data/From_Nottingham_with_Love/JRH_MotorCon_20100429_01_FORMARK.ds/MarkerFile.mrk 

tres=1/(D.fsample);
tim=1;
for tt=1:length(block_order),    
    x(tim:tim+block_length/tres-1,block_order(tt))=1;
    tim=tim+block_length/tres;
end;

figure;plot(D.time,x);

%%%%%%%%%
%% run glm to do regression against known finger tapping regressors

oat.source_recon.dirname=[workingdir '/subj1_results_beta_hmm' num2str(do_hmm) '_' oat.source_recon.method];
%oat.first_level.name=['wholebrain'];
oat.first_level.name=['first_level'];
oat=osl_load_oat(oat.source_recon.dirname, oat.first_level.name,'sub_level','group_level');

% GLM stuff:
oat.first_level.design_matrix=x';
oat.first_level.contrast{1}=[-1 0 1 0 0]'; % rest-left
oat.first_level.contrast{2}=[0 -1 1 0 0]'; % rest-right
oat.first_level.contrast{3}=[0  0 1 -1 0]'; % rest-both
oat.first_level.contrast_name{1}='rest-left';
oat.first_level.contrast_name{2}='rest-right';
oat.first_level.contrast_name{3}='rest-both';
oat.first_level.bc=[0 0 0];
oat.first_level.diagnostic_cons_to_do=1:3;
oat.first_level.cope_type='cope';
oat.first_level.tf_hilbert_freq_res=diff(oat.first_level.tf_freq_range);
oat.first_level.tf_method='hilbert';
oat.first_level.post_tf_downsample_factor=2; 
oat.first_level.tf_hilbert_freq_ranges=[13 30];
oat.first_level.time_moving_av_win_size=1;
oat.first_level.tf_hilbert_do_bandpass_for_single_freq=0;

oat.to_do=[0 1 0 0];

oat = osl_run_oat(oat);

%% output niis

S2=[];
S2.oat=oat;
S2.stats_fname=oat.first_level.results_fnames{1};
S2.first_level_contrasts=[1:3]; % list of first level contrasts to output
S2.stats_dir=[oat.source_recon.dirname '/' oat.first_level.name '_stats_dir'];
clear statsdir;
for ff=1:size(oat.first_level.tf_hilbert_freq_ranges,1),
    S2.freq_bin=ff;
    S2.stats_dir=[oat.source_recon.dirname '/' oat.first_level.name 'new_f' num2str(mean(oat.first_level.tf_hilbert_freq_ranges(ff,:))) '_stats_dir'];

    [statsdir{ff},times]=osl_save_nii_stats(S2);    
end;
% make sure you view results using fslview

contrast_num=3;
%runcmd(['fslview ' statsdir{1} '/tstat' num2str(contrast_num) '_2mm ' statsdir{2} '/tstat' num2str(contrast_num) '_2mm &']);
runcmd(['fslview ' statsdir{1} '/tstat' num2str(contrast_num) '_2mm &']);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Look at ROI time courses

%% First calculate an ROI mask 
oat.source_recon.dirname=[workingdir '/subj1_results_beta_hmm' num2str(do_hmm) '_' oat.source_recon.method];
oat.first_level.name=['first_level'];
oat=osl_load_oat(oat.source_recon.dirname, oat.first_level.name,'sub_level','group_level');
statsdir=[oat.source_recon.dirname '/' oat.first_level.name 'new_f' num2str(mean(oat.first_level.tf_hilbert_freq_ranges(ff,:))) '_stats_dir'];

% use FSL maths to threshold to create mask
con=3;
thresh=100;
runcmd(['fslmaths ' statsdir '/tstat' num2str(con) '_2mm -thr ' num2str(thresh) ' ' statsdir '/tstat' num2str(con) '_2mm_mask']);

% view the mask in fslview:
runcmd(['fslview ' statsdir '/tstat' num2str(con) '_2mm_mask &']);

%% Look at Hilbert Envelope timecourse averaged over ROI
oat.source_recon.dirname=[workingdir '/subj1_results_beta_hmm' num2str(do_hmm) '_' oat.source_recon.method];
oat.first_level.name=['first_level'];
oat=osl_load_oat(oat.source_recon.dirname, oat.first_level.name,'sub_level','group_level');

statsdir=[oat.source_recon.dirname '/' oat.first_level.name 'new_f' num2str(mean(oat.first_level.tf_hilbert_freq_ranges(ff,:))) '_stats_dir'];
oat.first_level.mask_fname=[statsdir '/tstat' num2str(con) '_2mm_mask'];
oat.first_level.design_matrix=x';
oat.first_level.contrast{1}=[-1 0 1 0 0]'; % rest-left
oat.first_level.contrast{2}=[0 -1 1 0 0]'; % rest-right
oat.first_level.contrast{3}=[0  0 1 -1 0]'; % rest-both
oat.first_level.bc=[0 0 0];
oat.first_level.diagnostic_cons_to_do=1:3;
oat.first_level.tf_hilbert_freq_res=diff(oat.first_level.tf_freq_range);

oat.first_level.tf_method='hilbert';
oat.first_level.post_movingaverage_downsample_factor=1; 
oat.first_level.post_tf_downsample_factor=2; 
oat.first_level.tf_hilbert_freq_ranges=[13 30];
oat.first_level.time_moving_av_win_size=1;
oat.first_level.tf_hilbert_do_bandpass_for_single_freq=0;
oat.first_level.do_glm_demean=0;
oat.first_level.name='roi';
oat.first_level.doGLM=0; % this turns the GLM off and just outputs the reconstructed data (that would be input into the GLM)

oat.to_do=[0 1 0 0];
oat = osl_run_oat(oat);

stats=osl_load_oat_results(oat,oat.first_level.results_fnames{1});
figure;plot(stats.glm_input_times,normalise(squeeze(mean(stats.glm_input_data,1))));

% compare to design matrix:
time_ind = ismembc2(stats.glm_input_times, D.time); %use undocumented built-in function
time_ind(time_ind==0) = [];

ho;plot(D.time(time_ind),x(time_ind,1),'r','LineWidth',2);
plot(D.time(time_ind),x(time_ind,2),'g','LineWidth',2);
plot(D.time(time_ind),x(time_ind,4),'k','LineWidth',2);
legend('beta','left','right','both');
plot4paper('time(secs)','power');

%%%%%%%%%
%% Look at ROI wideband time-frequency spectogram
% we first will need to redo beamformer with a wider band
oat.source_recon.dirname=[workingdir '/subj1_results_beta_hmm' num2str(do_hmm) '_' oat.source_recon.method];
oat.first_level.name=['first_level'];
oat=osl_load_oat(oat.source_recon.dirname, oat.first_level.name,'sub_level','group_level');

con=3;
oat.source_recon.mask_fname=[statsdir '/tstat' num2str(con) '_2mm_mask'];

% need new OAT dirname for new source recon
oat.source_recon.dirname=[workingdir '/subj1_results_wideband.oat'];

oat.source_recon.freq_range=[4 48];
oat.to_do=[1 0 0 0];
oat = osl_run_oat(oat);

%% now do first level for multiple freqs

oat.first_level.mask_fname=[statsdir '/tstat' num2str(con) '_2mm_mask'];
oat.first_level.tf_method='morlet';
oat.first_level.tf_morlet_factor=6;
if(isfield(oat.first_level,'tf_hilbert_freq_ranges') ),oat.first_level=rmfield(oat.first_level,'tf_hilbert_freq_ranges');end;
oat.first_level.tf_num_freqs=15;
oat.first_level.tf_freq_range=[4 48];
oat.first_level.tf_hilbert_freq_res=5;
oat.first_level.doGLM=0; % does not fit GLM and will output timeseries that would have been input into GLM
oat.first_level.post_tf_downsample_factor=2; 
oat.first_level.do_glm_demean=0;
oat.first_level.time_moving_av_win_size=10;

oat.to_do=[0 1 0 0];
oat = osl_run_oat(oat);

% load in results
stats=osl_load_oat_results(oat,oat.first_level.results_fnames{1});

figure;imagesc(stats.glm_input_times,stats.glm_input_frequencies,squeeze(mean(stats.glm_input_data,1))');colorbar;
axis xy;

% compare to design matrix:
ho;
x3=x(:,3);
x3(x3==0)=nan;
plot(D.time(time_ind),10*x3(time_ind),'k','LineWidth',5);ho;
legend('rest');
plot4paper('beta power','time(secs)');

